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INTRODUCTION 


The collision of tectonic plates involving large continental masses can produce significant crustal 
deformations, not only in the vicinity of the plate boundary, but also well into the interior of the 
plates. The best studied example of a contemporary continental collision is that involving the Indian 
and Eurasian plates, but other examples of both contemporary and ancient collisions are well known. 
The former category includes the collision between the Arabian and Eurasian Plates resulting in 
deformation in the Zagros Crush Zone, Anatolia, and elsewhere, and the incipient collision between 
Africa and the European portion of Eurasia. Ancient orogenic zones which display residual evidence 
of past continental collision include the Appalachian-Caledonian system, the Hercynian System in 
central Europe, and the Urals (Condie, 1982). Despite the attention given to both the general 
phenomena involved in continental collisions and to the geotectonics of specific collisions, the 
mechanical processes are still not well understood. In fact, controversy still exists about the correct 
kinematic description of the Indian-Eurasian collision. Nevertheless, the recent progress that has 
been made in understanding continental tectonics has sparked the emergence of several prototype 
models of collisional processes. The first of these models is the continental analog of the subduction 
of oceanic lithosphere under a continental plate near an oceanic trench. In the continental collision 
case crust from one plate overrides that of the other. The partial subduction of crust from one of the 
plates (which may have had a significant oceanic component as in the case of the Indian Plate) may 
be at a shallow angle with the possibility of considerable horizontal motion, or at a steep angle with 
implied penetration of the lithospheric plate into the mantle. In the case of the Indian-Eurasian colli- 
sion, there is little doubt that overthrusting of Asian continental crust onto lithospheric material from 
the oceanic portion of Indian plate occurred during the closing of the Tethys ocean. While localized 
thrusting occurs today in the vicinity of the Himalayas, it seems unlikely that this model can be used 
in the Tibet region, north of the Indus suture. One difficulty with the model is the fact that the 
buoyant Indian crust would tend to resist subduction. Both the absence of northly directed compres- 
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sional seismic activity, and the general absence of intermediate and deep depth seismic activity 
around Tibet argue against using this model to explain deformation within the Eurasian plate. 

A second model assumes a viscous deformation of one plate as a consequence of the ram-like 
penetration into this plate by a more rigid collisional partner. This model is not suitable for studying 
the thrusting responsible for the Himalaya orogen, but it has been used with considerable success by 
Tapponier and Molnar (1977) and others to explain the eastward extrusion of Tibet along the Altyn 
Tagh and Kun-Lun faults, the plateau-like thickening of Tibetan crust, and a number of other Asiatic 
tectonic features. 

A third model involving the shearing of thin layers of crust and their overthrusting one on top 
of another to accommodate the convergence between the plates has come to be known as “flake tec- 
tonics” (Oxburgh, 1972) and has been used in the study of both Appalachian and Alpian geological 
features. 

A model for Tibet in which the uniform topography is explained by having the plateau float like 
a hydrostatic head on a fluid has been discussed by several authors. Recently Zhao and Morgan 
(1985) expressed the view that relatively cold and rigid Indian crust is penetrating into warm, low 
viscosity Asiatic lower crust with Tibet rising in response to the induced hydrostatic pressure. These 
various models are not entirely mutually exclusive and many variants of each model can be 
advanced. 

In this paper we propose to explore in some depth various features of the punch/viscous sheet 
model of continental deformation. Our aim here is to understand the details of the model and to test 
how sensitive the model predictions are to the choices of numerical values for the parameters and to 
some of the implicit and explicit assumptions. We will limit the analysis to cases having simple 
geometries where our intuitive insight is sharper than might be the case for geologically truer condi- 
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tions. We refer the interested reader to the papers in the references for more detailed discussions of 
the model justifications and its applicability to specific contemporary environments. 


MODEL DESCRIPTION 

Figure 1 shows one version of the punch/deformable sheet model. By and large we will employ 
the mathematical development of this model that was derived by England and McKenzie (1982, 

1983). Original calculations based on this model used finite difference techniques which required that 
the punch boundary condition be expressed as a flux of material rather than a rigid indentation. 
Houseman and England (1985) improved the model by using finite element techniques with prescrib- 
ed boundary velocity conditions representing the punch motion. We use our own formulation of the 



Figure 1: Punch/viscous sheet model of continental collisions. The penetration of a rigid punch into 
a deformable sheet is represented by moving the boundary between points A and B north- 
ward at a uniform, constant velocity. Stress and deformation variables are averaged over 
the depth of the sheet. The sheet has a two layer density structure. The crustal density is 
Qc', this layer has a thickness, s. The mantle density is g m . The sheet is thin compared to 
its length and breadth, is subject to gravitational as well as collisional stresses, and 
deforms as an incompressible viscous fluid. 


3 



finite element approach with differences from the earlier work that will be introduced in the discus- 
sions to follow. The viscous sheet is assumed to be incompressible and in local isostatic equilibrium. 
The isostacy condition is maintained by Airy type compensation against a crustless ocean plate. The 
sheet thickness is assumed to be small compared to its length (the X or eastward direction) and 
breadth (the Y or northward direction). 


The deformation of the sheet is governed by a power law flow of the form: 

t . = BE("“') &y (!) 

where Ty and are components of the deviatoric stress and the strain rate tensor respectively, E is 
the second invariant of the strain rate tensor, and n is a power coefficient for linear (n = l) or 
nonlinear (n> 1) flow. The stress and strain rate variables are averaged over the thickness of the 
plate. The density structure of the plate is two layered. The crust has a density q c and a thickness, s 
(initial value s 0 ); the underlying mantle has a density p m . The initial thickness of the sheet is L; 
after collision its thickness at a point x,y is L + h(x,y) where h can be interpreted as a topographic 
height and we require h< <L everywhere. The isostatic condition relates the crustal thickness, 
topographic height, and densities: 

h = s(l -e c /e m ) ( 2 ) 


The depth averaged equilibrium equations provide the quasi-static equations of motion for the 
model (England and McKenzie; 1982, 1983), viz, 
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and 



where u and v are the velocities in the x and y directions respectively. The left hand sides of Equa- 
tions 3 and 4 contain dynamical quantities relating to the velocity field; while the right hand sides 
contain thickness and density variables which are quantities related to the gravitational stresses. The 
equilibrium conditions can be solved for the velocity field whenever the crustal thicknesses are 
known. In the general case, the continuity or conservation of mass condition provides a third equa- 
tion for the three unknown quantities: u(x,y,t), v(x,y,t), and s(x,y,t). For finite difference calcula- 
tions the continuity condition is best expressed in an Eulerian frame of reference with fixed positions 
for all nodal points. The Eulerian frame can also be used with finite element calculations, but a 
Langranian frame in which nodal points move with the material flow, is more convenient. In this 
frame of reference the continuity condition expressing conservation of crustal mass is: 

^ = — S (&xx ■*" &yy) = s &zz ^ 

The last form of the equality is a consequence of the incompressiblity condition, 

&XX +^yy + S zz = 0 (6) 

As we will now show the equilibrium equations have a useful and particularly simple form for a 
linear rheology. First we define a term, the gravitational strain rate, S g , which has strain rate units 
and can be used to express the contribution of gravity to the equilibrium conditions, 

p _ g ecO-ei/Cm) s 2 ( ? ) 
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The parameter B=2 -q where rj is the viscosity. The linear equilibrium equations can now be written 


or 


and 


or 


a d 2u . g2 u = _a_ [g e c (i - e c / e m )s 2 

3x 2 dy 2 3x3y 3x _ BL 


2 xx 
3x 


+- 


ag 


3x 


yy + 


ag 


ay 


xy = 


i ag 




2 ax 


4 d 2 v a 2 v _a 2 u_ = _a^ fg e c (i - g c / e m )s 2 

3y 2 3x 2 3y3x 3y _ BL 


2 



d& xx 

3y 


+ 


xy = 
5x 


1 

iV 


(8a) 


(8b) 


(9a) 


(9b) 


When the terms involving the gravitational strain are small, the equilibrium equations can be solved 
for the velocities; these velocities can then be used in the continuity equation. In the general case, 
the three equations must be solved simultaneously. Once the velocities and strain rate have been 
determined, deviatoric stresses can be calculated from the rheological law. To calculate total 
stresses, the thickness averaged pressure, P, must be known: 


p = ylg Gm L + Bg g ] - (t. 


xx T yy) 


( 10 ) 


Contraction is taken to be negative, extension positive, in our calculations. 

The parameters that enter explicitly into the problem are the rheological constants, B (or r; ) and 
n, the lithosphere thickness, L, and the densities, q c and g m . Initial conditions are imposed on the 
crustal thickness at the nodes, and boundary conditions are imposed on the velocities and (if desired) 
the thicknesses along the grid border. Constraints can be applied on the thickness of individual mesh 
points. The geometry of the plate is expressed through the set of grid parameters. 
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RESULTS 


A prototype case we have investigated involves a linear rheology with parameters chosen as in 
Table I. Our finite element grid consists of 361 nodes arranged in four node elements. The 324 
elements are initially squares with sides 200 km in length. The original grid is 3600 km in both 
length and breadth. A limited number of calculations were made with finer mesh grids to verify the 
general results presented here. We have taken configurations in which there is symmetry about the 
midpoint of the punch in Figure 1. The grid represents the right hand side of the deforming sheet; 
i.e., X=0 is the midpoint of the punch. The punch has a uniform velocity of 50 mm/yr in a direc- 
tion parallel to the y axis (northward) for values of X less than or equal to 1200 km. The bottom 
boundary of the sheet is held fixed for nodes at 1400 km and beyond. The top boundary is held 
fixed, but we have investigated the deformation for two lateral boundary conditions: one when the 
side boundary is held fixed and one where free slip is allowed (a stree free right hand boundary). 
Figure 2 shows velocity distributions at t=0 + , a time immediately after initial contact between the 
punch and sheet. In front of the punch, material flows generally northward to escape the indentation 

Table I: Numerical Parameters Used in Finite Element Model 

initial plate length: 3600 km 
plate breadth: 3600 km 
initial plate thickness, L = 100 km 
initial crustal thickness, s 0 = 35 km 
punch velocity: 50 mm/yr = 50 km/Myr northward 
punch width (center to edge): 1200 km 
crustal density, g c , = 2.80 gm/cm 3 
mantle density, g m , = 3.30 gm/cm 3 
B (=2 rj) = 1 x 10 24 Pa. sec 
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Further away, particularly to the northeast there is a tendency for material to flow toward the side of 
the grid; again in an attempt to move away from the zone of impact compression. As the side wall 
is approached the flow is outward for a free boundary, but must vanish for a fixed boundary. In the 
calculations shown on Figure 2, the side boundary is sufficiently far away that the velocities appear 
to vanish on the scale of the figure. However as later figures will show, the deformation at this 
boundary is not negligible. To the side of the indenter, a low pressure regime will develop as the 
punch progresses inward. As we will see in a moment, this becomes a region of crustal thinning and 
shearing or extensional movement. 


Some of the features of the collision model can be seen in more detail in Figure 3 which shows 
the distorted grid and crustal thicknesses after punch penetration has advanced one-third of the way 


S {Km) 

t = 24 Myr S 0 = 35 Km 



Figure 3: Grid distortions and crustal thicknesses (in km) after 24 million years. The grid elements 
are initially squares with sides of 200 km. The initial crustal thickness is 35 km. The right 
hand side is unconstrained. The element thickness values are the average of the values ob- 
tained at the four nodes of the element. 
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through the plate (at a velocity of 50 mm/yr, a distance of 1200 km is covered in 24 million years). 
The thickening of the crust is concentrated near the corner of the punch (but not at the edge due, in 
part, to the discrete nature of the finite element grid). The intense deformation in front of the corner 
of the punch is manifest by the distortion of the grid elements near that locale. Crustal thicknesses 
near the center of the punch are still substantial but not as large as those near the corner. As 
distance from the punch increases, there is a much more rapid decrease in crustal thickness for 
points located in front of the punch corner than for those in front of the punch center. To the side of 
the punch, the aforementioned low pressure zone manifests itself in the dilatation of the elements and 
the thinning of the crust. The bowing of the (free slip) right hand boundary results in crustal thin- 
ning in the boundary midsection and thickening toward the ends. Figure 3 should be compared to 
Figure 4 which shows the results when the lateral boundary is fixed. In general, crustal thicknesses 
are greater in the latter because material cannot expand beyond the original sides. In front of the 


S (Km) 

t .- 24 Myt S,, - 35 Km 



Figure 4: Same as Figure 3 except the right hand side is fixed. 
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punch, the qualitative features are very similar in the two cases. Therefore, when the aspect ratio of 
the punch diameter to sheet length is one-third or greater it appears that the deformation in front of 
the punch is controlled by the punch and not the lateral boundary conditions. Proceeding to the side, 
the differences between the two models become more pronounced as boundary effects become domi- 
nant. We have verified these conclusions by performing simulations with both smaller and larger 
punches. When the punch diameter exceeds fifty percent of the sheet length, the effects of the inden- 
tation are ubiquitous. 

The calculations shown in Figures 3 and 4 were terminated at 24 million years because the 
elements near the punch corner had become quite distorted. If the calculations were extended for a 
greater period, numerical instabilities would be introduced. We note, however, that the calculations 
could be extended in time by redefining the finite element grid and transforming all variables onto 
the new grid. Since the punch has already penetrated one-third of the way into the sheet by 24 
million years we chose not to introduce this complexity into our calculations. The choice of 24 
million years is arbitrary. By choosing larger distances between the nodal points, but maintaining all 
other parameters the integration could be extended further in time, but with some loss of spatial 
resolution. Spatial resolution could be improved by using a finer grid, but the calculations would re- 
quire greater computer time and storage. The calculations can also be extended in time by using a 
less abrupt edge for the punch as in Houseman and England (1985). However this alters the defor- 
mation pattern (see discussion below). 

Figures 5 and 6 focus attention on the instantaneous crustal thickening rates rather than the in- 
tegrated thickness over time. At t=0 + , the crustal thickening rates (which are shown as element 
values obtained by computing the mean of the four nodal values) have maximum values on the order 
of 1 mm/yr for the cases illustrated. The maximum rates are concentrated near the punch corner. 
After 24 million years the situation has changed somewhat. The region with moderate to large 
crustal thickening rates has grown larger and the location of the maximum has moved in front of the 
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Figure 5: Rates of crustal thickening at t=0+ for the same simulation as in Figure 3. The element 
rates are obtained by averaging the four nodal values for each element. 

punch. Concerning this latter observation we find that although gravity is acting on the topography 
to reduce the rate of uplift at the corner of the punch, this effect is not large. The extreme distortion 
of the elements around the punch corner make us cautious about speculating whether this effect is 
real or a limitation of the finite element procedure. 

The preceding discussion has emphasized crustal thickness and vertical movements. In fact, the 
original motivation for the punch/sheet model came from similarities between predicted and observed 
horizontal deformations. Figures 7 and 8 show computed values of the horizontal components of the 
strain rate tensor. 

The analogy between some of the features observed from this model and those observed in Asia 
was first discussed by Molnar and Tapponier (1975; see also Tapponier, et al., 1982). Among the 
most important characteristics are the crustal thickening in front of the punch, the shearing and 
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Figure 6: Same as Figure 5, but t= 24 million years. 

lateral extrusion of material to the east beyond the corner of the punch, the termination of this shear- 
ing by either compression (in the case of fixed boundaries) or extension (free boundaries), and shear 
to the east of the punch. These phenomena are all observed to some degree in Asia, and seem to be 
a consequence of the collision with India. 


Since north-south compression dominates the deformation in front of the punch it is interesting 
to take a closer look at the behavior of this quantity. The compression rate is greatest at the corner 
of the punch, and decreases with distance away from the punch. As shown in Figure 9, however, 
the compression in front of the center of the punch increases out to a distance of about 1000 km 
from the impact line, then decreases. The position of this extremum shows little variation with time 
if the origin is moved northward along with the advancing punch. A similar qualitative behavior can 
be deduced from the analytical expressions of Houseman, et al., but the quantitative behavior seen 
here is quite different. The analytical equations predict that the strain rate maximum should occur at 
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DISTANCE NORMAL TO PUNCH (10 2 Km) 

Figure 9: North-south strain rate versus distance from the punch. 

x = lambda/(6 pi) where lambda is the wavelength of the (sinusoidal approximation to the) punch. 
Choosing lamdba to be about 5000 km corresponding to the velocity dropping from its peak value to 
zero over a 1/4 wavelength interval of 1250 km gives a predicted peak at 265 km from the punch. 
This analytically predicted peak occurs much closer to the punch than is observed in the numerical 
simulations. Apparently significant inaccuracies are introduced in the analytical calculations by the 
approximations that are employed. In particular, the assumption of a sinusoidal variation in the 
velocity distribution on the x axis produces a region of southward as well as northward flow. The 
transverse components of stress and flow are also ignored in the analytical treatment. The near in- 
variance of the location of the strain rate extremum appears to be a feature of a linear rheology. 
Preliminary calculations using a power law rheology with n=3 indicate that the extremum migrates 
away from the punch much more rapidly in the nonlinear case. A related modeling result by Melosh 
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PRINCIPAL HORIZONTAL STRAIN RATES (t = 0 + ) 



Figure 10: Principal strain rates at t=0+ for free lateral boundary problem. 

(1976) shows that postseismic stress propagates more rapidly into the interior of a layered nonlinear 
(as opposed to linear) viscoelastic earth. 

An example of the principal strain rates that are deduced from the model are shown in Figure 
10. The modes of deformation in different locales can be classified by applying a generalization of 
the Anderson brittle failure criteria to the viscous flow occurring here. Accordingly a compressive 
(thrusting) regime dominates where the vertical strain rate is the largest of the principle strain rates, 
a shearing regime exists when the vertical strain rate is the intermediate value of the three principal 
values, and an extensional (normal) regime dominates when the vertical strain is the least of the 
principle values. The classification is shown for two different times in Figure 11. The portion of the 
sheet dominated by compression increases with time and the shearing zones decrease in areal extent. 
In particular, the zone of left lateral shearing occurring to the immediate northwest of the punch is 
displaced northward as the punch progresses into the sheet. Of course, the classification of different 
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Figure 11: Strain rate regimes classified according to dominant deformation mechanism. Free lateral boundary 

a. t=0 + . 

b. t=24 Myr. 



regions is strongly influenced by boundary conditions. Not only is the extensional behavior on the 
right side dependent on the stress free boundary conditions there, but also it would be reasonable to 
expect that shearing would be more important near the northern boundary if roller-type boundary 
conditions permitting eastward sliding were imposed. 

Another integrated quantity which naturally emerges from the numerical calculations is the 
change in the areal extent of the elements adopted for the finite element grid. For the problem with 
a free lateral boundary, areal contraction occurs to the north and east of the punch whereas areal ex- 
pansion occurs to the southwest and extreme east as shown in Figure 12. For the problem with a 
fixed lateral boundary, the region of contraction is larger (contraction occurs at the lateral boun- 
daries) and that of expansion, smaller. 

The preceding discussion has focused attention on the general behavior of the punch/viscous 
sheet model and compared results using different boundary conditions. Although an examination of 
all reasonable variations in model parameters would require considerably more space than can be 
allotted to this paper, some insight into the sensitivity of model results to parameter variations can 
be presented here. Viscosity, crustal densities, and lithospheric thickness enter the calculations 
through the right hand side of equations 3, 4, 8 and 9. In most of our calculations with a linear 
rheology we have found that this term has little importance. In fact most of the observed horizontal 
deformation patterns remain substantially unaltered if the crustal and mantle densities are set equal. 
For the densities shown in Table I, the gravitational term tends to be small for B values (= two 
times viscosity) greater than about 1 x 10 23 23 Pa. sec unless the lithosphere is very thin or gra- 
dients in crustal thickness are very large. For example, the grid distortions and crustal thicknesses 
for B = 1 x 10 22 22 Pa. sec are shown in Figure 13. The reduced crustal thicknesses are due to the 
inability of the sheet to maintain significant vertical straining against gravitationally induced stresses. 

The distribution of velocities on the punch has a profound effect on the results. Although this 
statement is intuitively obvious, it has considerable importance for comparisons between the present 
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Figure 13: Low viscosity result analagous to Figure 3. B = 1 x 10 22 Pa. sec, or 1/1 00th that of 
the earlier results. x 

calculations with the earlier finite difference work of England and McKenzie (1982, 1983). In the 
earlier work, the punch could not be represented by a moving boundary, so a flux boundary condi- 
tion was imposed on the bottom boundary of the grid. The calculations were performed on a fixed 
grid in an Eulerian frame of reference. As a consequence the grid of nodal points was fixed in space 
rather than moving with the deformation. The flux boundary conditions was a flow of material at a 
velocity of 50 mm/yr northward between the boundary center and a point half-way across the punch. 
The velocity distribution then decreased to zero at the punch edge via a cosine squared law. In the 
early stages of the collisions, the vertical straining, compression, and crustal thickening are concen- 
trated at the edge of the punch. The observed features are qualitatively similar to those seen with the 
finite element, moving boundary calculations except, of course for the deformation observed south of 
the boundary in the flux calculations. As time advances, a punch boundary can be defined by 
calculating the position of material that was initially on the grid boundary when the collision began. 
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When this is done, as in Figure 14, one finds that the punch becomes progressively more sheared as 
it deforms according to the same rheological law as the sheet. This deformation includes a rounding 
of the velocity distribution. The center of the boundary becomes the leading point and the sharp cor- 
ner at the punch edge is eliminated. As a consequence crustal thicknesses become concentrated near 
the middle of the leading edge of the advancing material . The pattern of horizontal straining is also 
altered to some degree. These results suggest several important features about the deformation. The 
concentration of the most intense deformation about the punch corner requires the maintenance of a 
sharp corner throughout the duration of the collision. If the punch is to represent the continental por- 
tion of a tectonic plate, then the continental material must remain fairly rigid and unsheared. If, 
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Figure 14: Crustal thicknesses and punch boundary obtained from flux boundary condition problem, 
t =32 Myr. A flow northward of 50 mm/yr occurs between X=0 and X = 800 km; the 
velocity then decreases to 0 at X = 1600 km according to a cosine squared law. The finite 
element grid is fixed and a Eulerian. analysis is used. The solid line shows the position of 
material that was on the southern grid boundary at t=0. This punch boundary has pro- 
gressively sheared both due to the velocity distribution and the non-rigid character of the 
punch. Crustal thickness has a maximum at X=0. 
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however, the continental ram does undergo significant shear, then the deformation within the sheet 
becomes concentrated in front of the leading point of the punch. 

SUMMARY 

Our modeling efforts are based almost exclusively on the geophysical and mathematical founda- 
tion provided by the work of England, McKenzie, Houseman, and Sonder. Therefore, some of the 
results presented here for the fixed lateral boundary conditions will qualitatively agree with the finite 
element results of these workers. The new features introduced here include a free lateral boundary, a 
rigid punch (as opposed to a tapered and sheared punch used in earlier work), and a more detailed 
analysis of several aspects of the model. Particularly important features found in our calculations are 
the sensitivity of compression to position from the punch, the relative absence of strain propagation 
away from the punch for a linear viscous flow (as opposed to the propagation with a power law 
flow), the termination of shearing by compression or extension depending on boundary conditions, 
and the effect of punch rigidity on deformation near the corners and the center of the punch. New 
results have been presented on crustal thickening rates (which approach 1 mm/yr for the parameters 
used herein), changes in the area of the mass elements, and the classification of deformation 
regimes. 
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